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Abstract. We analyze the dynamics of a 3D granular packing composed of particles of irregular polyhedral 
shape confined inside a rectangular box with a retaining wall subjected to horizontal harmonic forcing. The 
simulations are performed by means of the contact dynamics method for a broad set of loading parameters. 
We explore the vibrational dynamics of the packing, the evolution of solid fraction and the scaling of dy- 
namics with the loading parameters. We show that the motion of the retaining wall is strongly anharmonic 
as a result of jamming and grain rearrangements. It is found that the mean particle displacement scales 
with inverse square of frequency, the inverse of the force amplitude and the square of gravity. The short- 
time compaction rate grows in proportion to frequency up to a characteristic frequency, corresponding to 
collective particle rearrangements between equilibrium states, and then it declines in inverse proportion to 
frequency. 

PACS. 83.80.Fg Granular solids - 45.70.Cc Static sandpiles; granular compaction 



1 Introduction 

The dynamics of dense granular materials subjected to 
vibrations involves collective phenomena resulting from 
kinematic constraints (steric exclusions, weak spatial di- 
mensions, . . . ) and energy dissipation pQ. Well-known ex- 
amples of the vibration- induced phenomena are compaction, 
convective flow, size segregation and standing wave pat- 
terns at the free surface gHUMSUZUH] • Three different 
states can be distinguished depending on the intensity and 
frequency of vibrations: 1) Gas-like or fluidized state: The 
rate of energy input is such that there are no enduring 
contacts between particles and the material behaves as a 
dissipative gas [ MTT1I17] . 2) Solid-like state: Vibrational 
energy propagates through the network of enduring con- 
tacts between particles and the material undergoes slow 
rearrangements and progressive compaction 12, 2Tl llO[[2"0] . 
3) Liquid-like state: Both particle migration and enduring 
contact networks are involved in the dynamics and various 
collective effects can be observed [3l [T8lH3l [5lfT9] . 

Vibro-compaction is the main feature of the solid-like 
state. Most investigated systems are unconfined granular 
beds (with a free surface) subjected to vertical vibrations. 
The vibrations behave as a source of randomness allow- 
ing the system to explore metastable configurations and 
to reduce its potential energy. There is, however, another 
mechanism which prevails in the case of confined granu- 
lar materials. It is well-known that under cyclic straining, 



a granular material accumulates plastic deformation and 
the solid fraction tends to a maximum value depending on 
the material. This phenomenon is sometimes called "gran- 
ular ratcheting" due to the irreversible character of com- 
paction under cyclic loading [IIIITIjIITo] . In both confined 
and unconfined geometries, the solid fraction evolves as 
a logarithmic function of the number of cycles. In most 
work reported on vibrated granular media, the collective 
dynamics of the particles and the influence of various pa- 
rameters related to the material or the driving system have 
not been investigated in all details. Moreover, in nearly 
all studies, spherical or nearly spherical particles in 3D or 
disks or polygons in 2D have been used [14ll331 [T6] . 
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In this paper, we present a numerical investigation of 
the dynamics and short-time compaction of a system of ir- 
regular polyhedral particles confined inside a rectangular 
box with a retaining wall subjected to horizontal harmonic 
loading. This system is different from nearly all experi- 
mental systems investigated under horizontal vibrations 
since the packing remains confined inside a box, so that 
the gravity plays little role during the inward motion of 
the retaining wall. This geometry is similar to that used in 
various industrial applications such as the casting of fresh 
concrete where efficient vibro-compaction of dry and wet 
granular materials represents a crucial issue [j22,123|. The 
tamping operation on railway ballast is another interest- 
ing case where vibrating tamping bars are used to restore 
the initial geometry of the track distorted as a result of 
ballast settlement [2"4fl25]. With the increase of commercial 
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speed, a better understanding of the physics of compaction 
is important for long time stability of ballast. 

We are interested here in the evolution of the packing 
in the course of harmonic loading, the short-time com- 
paction (during the first cycles) and the scaling of the dy- 
namics with loading parameters. We used discrete-element 
numerical simulations by means of the contact dynam- 
ics method in 3D with rigid irregular polyhedral particles 
p6l[27l[28ll34ll29] . The system is explored for a broad set 
of loading parameters including the frequency and ampli- 
tude of the harmonic driving force. A similar study was re- 
cently performed in 2D with polygonal particles L 33] . The 
passage from 2D to 3D from a purely numerical point of 
view involves numerical handling of particles of polyhedral 
shape and a higher numerical efficiency making 3D simu- 
lations over many cycles and for varying parameters pos- 
sible. In this paper, we will revisit the same phenomenol- 
ogy as in 2D for irregular polyhedral particles. This shows 
that in granular materials the dynamics is not sensitive 
to space dimensionality although the influence of particle 
shape and the details of structural organization can only 
be appreciated in a 3D geometry. 

We begin with an introduction to 3D contact dynamics 
method as applied to polyhedral particle shapes and the 
numerical procedures. Then, we present in three sections 
the dynamics of the packing, the evolution of solid fraction 
and scaling with loading parameters. 

2 Numerical method 

The simulations were carried out by means of the contact 
dynamics (CD) method with irregular polyhedral parti- 
cles. In this section we present the properties of this nu- 
merical method and compare it to a more classical numeri- 
cal approach molecular dynamics (MD). Then we presents 
the principle of the contact detection and numerical pa- 
rameters which have been used for this study. 

2.1 Contact dynamics 

The CD method is based on implicit time integration and 
nonsmooth formulation of mutual exclusion and dry fric- 
tion between particles in case of contact 26 , 27, 28 , 34j . The 
equations of motion for each particle are formulated as 
differential inclusions in which velocity jumps replace ac- 
celerations [37] . The unilateral contact interactions and 
Coulomb friction law are represented as set-valued force 
laws according to convex analysis. The implementation of 
the time-stepping scheme requires the geometrical descrip- 
tion of each potential contact in terms of contact position 
and its unit normal vector. 

At a given step of evolution, all kinematic constraints 
implied by enduring contacts and the possible rolling of 
some particles over others are simultaneously taken into 
account, together with the equations of dynamics, in order 
to determine all velocities and contact forces in the sys- 
tem. This problem is solved by an iterative process (non- 
linear Gauss-Seidel method) which consists of solving a 



single contact problem, with other contact forces being 
treated as known, and iteratively updating the forces un- 
til a given convergence criterion is fulfilled. The method 
is thus able to deal properly with the nonlocal character 
of the momemtum transfers — resulting from the perfect 
rigidity of particles in contact. The CD method makes no 
difference between smooth evolution of a system of rigid 
particles during one time step and nonsmooth evolutions 
in time due to collisions or dry friction effects. 

The MD-like methods are based on regularization schemes 
where impenetrability is approximated by a steep repul- 
sive potential and Coulomb's law by a viscous- or elastic- 
regularized friction law, to which smooth computation 
methods can be applied. In this case, the choice of vis- 
cous or elastic parameters has to be specified depending 
on the particle shapes. This regularization implies small 
time steps in order to ensure numerical stability whereas 
the implicit time integration method, inherent in the CD 
method, is unconditionally stable. The uniqueness is not 
guaranteed by CD approach for perfectly rigid particles in 
absolute terms. However, by initializing each step of calcu- 
lation with the forces calculated in the preceding step, the 
set of admissible solutions shrinks to fluctuations which 
are basically below the numerical resolution. In MD-based 
simulations, this "force history" is encoded by construc- 
tion in the particle positions. 

For our simulations, we used the LMGC90 which is 
a multipurpose software developed in Montpellier, capa- 
ble of modeling a collection of deformable or undeformable 
particles of various shapes (spherical, polyhedral, or polyg- 
onal) by different algorithms 34J. 

2.2 Simulation of polyhedral particles 

The treatment of a contact interaction between two par- 
ticles requires the identification of the contact zone and a 
plane (in 3D), the so-called "common plane". Obviously, 
for rigid particles it is possible to define this contact zone 
by a finite set of points. Before applying the contact de- 
tection algorithm between a pair of particles of irregu- 
lar shapes (see Fig. [5]), a "bounding box" method is used 
to compute a list of particle pairs potentially in contact. 
Then, for each pair, the first step is to determine if an over- 
lap exists through a 3D extension of the "shadow overlap 
method" [241134] , Several algorithms exist for overlap de- 
termination between convex polyhedra 35,34,36 . In the 
case of an overlap, the contact plane is determined by 
computing the intersection between the two particles. 

The contacts between polyhedral particles belong to 
different categories, namely face-face, edge-face, vertex- 
face, edge-edge, vertex- vertex, vertex-edge; see Fig.[TJ The 
vertex-vertex and vertex-edge contacts are rare, i.e. in 
practice of "zero measure" . For each case we determine 
one, two or three contact points which provide a good 
description of the contact zone. In this paper, the vertex- 
edge and edge-edge contacts are referred to as "simple" 
contacts. On the other hand, the edge-face and face-face 
contacts are treated as "double" and "triple" contacts, re- 
spectively, as their representation involves 2 and 3 distinct 
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Fig. 1. Different types of contacts between two polyhedra. 



points on the common plane. The detection procedure is 
fairly rapid and allows us to simulate large samples com- 
posed of polyhedral particles. 



2.3 Numerical samples 



Our numerical samples are composed of rigid polyhedral 
particles with shapes and sizes that represent those of bal- 
last grains (Fig. [2]). Each particle has at most 70 faces and 
37 vertices and at least 12 faces and 8 vertices. A sample 
contains nearly 1200 particles. The particle size is char- 
acterized as the largest distance between the barycenter 
and the vertices of the particle, to which we will refer as 
"diameter" below. We used the following size distribution: 
50% of diameter d min — 2.5 cm, 34% of diameter 3.75 cm, 
16% of diameter d max = 5 cm. The particles are initially 
placed in a rectangular box and compressed by downward 
motion of the upper wall at zero gravity; see Fig. [3l Then, 
the gravity is set to g and the upper wall is raised 1 cm 
and fixed. The coefficient of friction between the particles 
and with the horizontal walls was fixed to 0.4, but it was 
at the vertical walls. The coefficient of restitution between 
particles was fixed to zero because of the high solid frac- 
tion of the samples. One of the walls is allowed to move 
horizontally (x direction in Fig. [3]) and subjected to a har- 
monic driving force. All other walls are immobile. For all 
simulations the time step was 2.10 _4 s. 




Fig. 2. Examples of polyhedral shapes used in the simulations. 




Fig. 3. A snapshot of the packing inside a box with a free 
wall over which the driving force f(t) is applied along the x 
direction. 



3 Vibrational dynamics 

The free wall is subjected to a harmonic force as a function 
of time, 



fit) 



if max ~t~ fjnin) (frr 



fmin ) . , / , 
SlUOjt, (1) 



where f max and f m i n are the largest and lowest compres- 
sive (positive) forces acting on the wall. The first term rep- 
resents the mean confining force modulated by the second 
term. At t = 0, the external force / = (f max + /min)/2 
causes the inward motion (contraction) of the free retain- 
ing wall. Jamming occurs when the gap left between the 
upper wall and the free surface of the packing is filled. If 
fmin is above the (gravitational) force exerted by the par- 
ticles on the free wall, / will be large enough to prevent the 
wall from backward motion (extension) during the whole 
cycle. Then, the granular material is in "passive state" 
and the major principal stress direction is horizontal |37J. 
On the other hand, if fmax is below the force exerted by 
the particles, / will never be large enough to prevent the 
extension of the packing. This corresponds to the "active 
state" where the major principal stress direction remains 
vertical. 
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Fig. 4. The evolution of the displacement x of the free wall 
(up) in response to harmonic loading (down). 
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Fig. 5. A zoom on the displacement of the free wall as a func- 
tion of time for a single period (up) in response to harmonic 
force (down). 



In all other cases, both contraction and extension occur 
during each period, and the displacement Ax of the free 
wall will be controlled by f m in an d fmax- Without loss 
of generality, we set f min = 0. This ensures the largest 
possible displacement of the wall in the active state. Four 
different values of f max were tested, ranging from 2.10 3 N 
to 10 4 N. 

We first consider the trajectory x(t) of the free wall 
which reflects the dynamics of the particles in the box 
in response to harmonic forcing. Figure 2] shows x(t) for 
frequency v = 5 Hz over a time interval At — 1 s. We can 
observe a fast initial contraction (t < 0.1 s) followed by 
slow contraction over four periods. The initial contraction 
is a consequence of the gap left between the free surface of 
the packing and the upper wall. The subsequent periodic 
motion of the wall takes place around this jammed state 
and will be studied below. 




Fig. 6. Instantaneous particle velocity field in the passive 
state, i.e. during inward motion of the free wall. 
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Fig. 7. Force f g exerted by the particles on the free wall as a 
function of displacement x. 



A zoom on a single period is shown in Fig. [5] It begins 
at the jamming position x — Xj corresponding to the jam- 
ming position reached at the end of the preceding period. 
The motion of the wall begins (point a in Fig. [5]) only 
when the applied force / declines near to its minimum 
fmin = 0. The maximum displacement Ax max occurs at a 
later time 6t r (point b). From a to b, the force exerted by 
the packing on the free wall is above the applied force, so 
that the wall moves backward (extension). In this phase, 
the packing is in an active state. The inverse situation pre- 
vails from b to c where the particles are pushed towards 
the box (contraction). Then, the packing is in a passive 
state. The new jamming position x'j is below the jamming 
position Xj reached at the end of the preceding period. 
The difference Xj — x'j represents the net compaction of 
the packing over one period. The particle velocity field is 
not a simple oscillation around an average position. The 
particles undergo a clockwise convective motion in the cell 
as shown in Fig. [6] 
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Fig. 8. Normal forces in the passive (a) and active (b) states in 
a section of the packing parallel to the xz plane. The segments 
connect particle centers with a thickness proportional to the 
normal force. 



Figure [7] shows the horizontal force f g exerted by the 
packing on the wall as a function of x over four periods. In 
the active phase, f g grows slightly with x. In the passive 
phase, it grows faster and almost linearly as x decreases. 
The vertical line corresponds to the jammed state where 



f g decreases with f at x 



We also observe two tran- 



sients : 1) unjamming and the onset of the active state, 2) 
jamming from the passive state. Inside the packing, the 
contact forces evolve between a fully jammed state, where 
horizontal force chains dominate (Fig. [HJa)), and the ac- 
tive state, where vertical gravity-induced chains can be 
observed (Fig. Mb)). 



4 Granular ratcheting 

In our system, the solid fraction p of the packing increases 
due to horizontal vibrations. This accumulation of plas- 
tic strain under oscillatory loading is sometimes called 
"granular ratcheting" [14,15]. To evaluate p, we consider 
a control volume inside the box. The initial value of the 
solid fraction is 0.50. Figure [HI shows the evolution of the 
variation Ap of solid fraction for several periods. An ini- 
tial compaction of 2% is followed by oscillations with a 
small increase of Ap in each period. The initial compaction 
should be attributed to the initial state of the packing is 
not yet fully confined. We use po = 0.51, reached after a 




0.4 0.6 

t(s) 

Fig. 9. Evolution of the solid fraction Ap from the initial state 
as a function of time over several periods. 



time laps of 0.2 s, as the reference value for the evolution 
of solid fraction. The relative compaction of the packing 
is given by Ap/ p^. The compaction rate fj over several 
periods and for a total time interval At is 



■ - 

V = poAt' 



(2) 



The compaction of the packing slows down logarithmically 
at long times [18]. But, the short-time compaction can 
well be approximated by a linear function with a constant 
compaction per period Ap\. Then, we have 



n 



Api 
Po 



(3) 



For v = 5 Hz and f max = 6 10 3 N, we have f] ~ 0.025 s _1 . 
This rate is faster in 3D compared to 2D simulations for 
the same frequency [33]. 

It is important to note that compaction occurs in the 
active state, i.e. during the extension of the packing. This 
is shown in Fig. [10] where the variation Ap of the solid 
fraction is plotted as a function of x. The solid fraction 
increases during extension (increasing x) and decreases 
during contraction (decreasing x). 

Granular ratcheting has been investigated by numeri- 
cal simulations showing that the anisotropy of critical con- 
tacts, where the friction force is fully mobilized, plays an 
important role 14.15J. Quasi-static cyclic shearing also 
leads to cumulative compaction of a granular material at 
low strain amplitudes [351l38| . At large amplitudes, the 
compaction is followed by decompaction (dilation) and no 
net compaction can be observed over a full cycle. In our 
system, compaction is a consequence of unjamming and 
it is pursued during the whole active state. Decompaction 
takes place in the passive state, but it is cut short by fast 
jamming. The outcome of a full cycle is thus a net com- 
paction of the packing. 
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Fig. 10. Variation Ap of the solid fraction from the initial 
state as a function of the displacement x of the free wall. 

5 Influence of loading parameters 

We performed a series of simulations for frequencies v 
ranging from 1 Hz to 60 Hz and for a total time of 1 
s. All simulations yield similar results both for dynamics 
and compaction. Moreover, a simple dimensional analysis 
leads to the collapse of the data on a single plot. Indeed, 
the frequency sets the time scale r = v~ x . Force scales are 
set by the largest driving force f max in the passive state 
and the particle weights mg as well as the smallest driv- 
ing force fmin in the active state. Hence, dimensionally, 
for fixed values of mg, f m in and f ma x, all displacements 
are expected to scale with v~ 2 and all velocities with v^ 1 . 

This scaling is shown in Fig. [TT] where the phase space 
trajectory is shown for v = 5 Hz and v = 10 Hz without 
scaling and after scaling the displacements Ax by v~ 2 
and the velocities v by v~ x . We see that the data from 
both simulations collapse nicely on the same trajectory 
after scaling. Figure [12] shows the maximum displacement 
Ax max in the active state and the maximum velocity v max 
in the passive state as a function of v. The corresponding 
fits by v~ 2 and v^ 1 are excellent. 

The role of force parameters rag, f m in and f max is less 
evident. Since we have / m , n = 0, we expect Ax max to be 
dependent on the ratio mg/ f max representing the relative 
importance of the gravitational to driving forces. Indeed, 
our data show that Ax max varies as J^; Fig. [13] On 
the other hand, the mass ratio m w /m, where m w and m 
are the mass of the free wall and the total mass of the 
particles, must control the inertia and thus the maximum 
displacement of the wall. Our simulations with different 
values of m w show that Ax max varies as mj(m + m w ). 

Hence, we propose the following expression for the scal- 
ing of displacements with loading parameters: 



C 



mg 

fmax 



(4) 



where C is a dimensionless prefactor. Fig. ll5l shows Ax rnax 
as a function of (mg) 2 / [(m + m w )(f rnax i/ 2 )] from different 
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Fig. 11. Phase space trajectories for two frequencies without 
scaling (a) and with scaling (b) of the displacements and ve- 
locities with respect to the frequency. 



simulations with different values of v, f max , g and m w . 
We see that the data are in excellent agreement with Eq. 
[U The prefactor is C ~ 0.01. This scaling is the same 
as in 2D simulations with a material constant C ~ 0.05 
for polygonal particles [33] . Let us also remark that Eq. [4] 
predicts that Ax max varies as g 2 . This prediction agrees 
well with our simulation data shown in Fig. fT?] for four 
different values of g. 



6 Compaction rates 

We now come back to granular ratcheting and we would 
like to evaluate short-time compaction rates as a function 
of frequency. According to Eq. [3] the compaction rate fj 
varies linearly with the frequency v if the total compaction 
per period Ap\ is independent of v. Fig. [1^] shows f\ as a 
function of v. We see that only at low frequencies, rj in- 
creases linearly with v. Beyond a characteristic frequency 
i> c , 77 declines with v. The largest compaction rate i] max oc- 



curs for v 



10 Hz. The corresponding time r c 



Emilien Azema et al.: Short-time dynamics of a packing of polyhedral grains under horizontal vibrations 



7 



0.002 



§ 0.001 



3 



0.000, 
0.08 

0.06 

,0.04 

0.02 

0.00 



• data 



10 20 30 40 50 60 

v(Hz) 



II 1 1 


1 1 






• data 

-l 

~v 






(b) 




> 












1 1 1 


• • . 

i i 





10 20 30 40 50 60 

v (Hz) 

Fig. 12. Maximum displacement Ax ma x (a) and the maximum 
velocity v max (b) as a function of frequency v. 
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Fig. 13. Scaling of the maximum displacement Ax ma x with 
the force amplitude fmax- 
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Fig. 14. Scaling of the maximum displacement Ax ma x with 
gravity g. 
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Fig. 15. Scaling of the maximum displacement Ax max with 
loading parameters from simulations with different values of 
the frequency v (triangles), the force amplitude f m ax (circles), 
the gravity g (squares), and for the mass m m (diamonds) of 
the free wall. The inset shows the plot near the origin. 



represents the characteristic time for the relaxation of the 
packing. In the active state, the packing needs a finite rear- 
rangement time to achieve a higher level of solid fraction. 
As long as the vibration period r = v~ x is longer than t c , 
the packing has enough time to relax fully to a more com- 
pact equilibrium state. But, when the period t is below r c , 
the relaxation is incomplete so that Ap\ < Ap max , where 
Ap max is the largest compaction between two periods. 

It is expected that Ap\ should follow the same scaling 
with the frequency as the displacement of the retaining 
wall, i.e. Api cx Ap max v~ 2 . This is because the volume 
change AV is proportional to Ax. Hence, from Eq. [3] and 
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Fig. 16. The compaction rate r\ as a function of the frequency 
(circles) fitted by two different functions; see text. 



imposing the continuity at v — v c , we get 



Apr, 



Pa 



V V < v c 



(5) 



Fig. [TBI shows this prediction together with the data points. 
We see that, although v c is the only fitting parameter, the 
compaction rate rj is well adjusted by Eq.[5l The prefactor 
Ap max /po is ~ 0.005, corresponding to Ap max ~ 0.0025. 
In spite of the sharp transition at v = v Cl it is convenient 
to construct a single expression containing the correct be- 
havior both at low and high frequencies. As in 2D for 
polygon packings [33], the following fitting form provides 
a good approximation as shown also in Fig. [16] (fitting 
form 2): 

V - 



1 



(0) 



The characteristic time r c = 0.1 s is of the same or- 
der of magnitude as the time required for one particle 
to fall down a distance equal to its diameter. Obviously, 
the above findings concern only short-time compaction 
(At < 1 s). At longer times, fj declines with time, but 
the scaling with frequency according to Eq. [5] is expected 
to hold at each instant of evolution of the packing. 



7 Conclusion 

In this paper, the contact dynamics method was employed 
to simulate and analyze the dynamics of a system of poly- 
hedral particles subjected to horizontal harmonic forcing 
of a retaining wall. Our system is devoid of elastic ele- 
ments and, hence, the behavior is fully governed by parti- 
cle rearrangements. Moreover, it involves a jammed state 
separating passive (loading) and active (unloading) states. 
Dimensional analysis was used to scale the displacements 
with the frequency of oscillations. It was shown that the 



data collapse by scaling the displacements by the inverse 
square of frequency. We also studied the scaling with con- 
fining force and particle weights. 

Granular ratcheting under horizontal vibrations was 
investigated. During each vibration period a small com- 
paction of the system occurs during unloading, i.e. upon 
sample extension, followed by decompaction upon con- 
traction. The compaction rate increases linearly with fre- 
quency up to a characteristic frequency and then it de- 
clines in inverse proportion to frequency. The characteris- 
tic frequency was interpreted in terms of relaxation time 
of the packing under its own weight during the unloading 
phase. 

The similarity of the phenomenology of vibrational dy- 
namics and compaction at short times in the 3D system 
of polyhedral particles with that of a 2D system of polyg- 
onal particles suggests that space dimensionality plays a 
minor role in granular dynamics. The characteristic times 
and compaction rates are slightly different but the scal- 
ing behavior and the functional dependence of the com- 
paction rate with frequency are the same. A comparison 
with spherical particles would be interesting in order to 
highlight the effect of particle shapes on these parameters. 
On the other hand, the characteristic time appears as a 
crucial parameter for the compaction rate and it merits 
further investigation as a function of various control pa- 
rameters of the system. 

The authors would like to thank specially F. Dubois for 
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